plot density of LAI values across years

this is for needleleaf - canary island pine

## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

I like to start with these density plots - really showing the distribution of outcomes for the different parameter sets. I started with this small value because it runs faster to compile these plots and get an idea of if there are any differences. This plot tells me there are a solid set of runs that included reasonable LAI for all three years, but not that there was recovery by 2017.

To compare how much each parameter set changes from 2011-2017, I plotted those against each other. Almost every run has a greater value in 2011.

I also took the running mean, for every 5 years, and plotted that for the entire output (2005-2025), shown below. There is a general downward trend for both LAI and plant C from the start.

plot running mean, and filter out based on LAI values in 2011

data =  read.csv("~/Google Drive/My Drive/patches/R/data/weighted_tree_summary.csv")
quag <- data[data$class == 34,]
low_lai <- quag$mean_wt_lai - 2*quag$sd_wt_lai
high_lai <- quag$mean_wt_lai + 2*quag$sd_wt_lai
low_carb <- quag$mean_wt_carb - 2*quag$sd_wt_carb
high_carb <- quag$mean_wt_carb + 2*quag$sd_wt_carb

test2 <- read %>% group_by(year, run) %>% 
  summarize_at(vars(plantc, lai), list(mean)) %>%
  arrange(run, year) %>%
  mutate(ma2=rollapply(lai, 5, mean, fill=NA)) %>%
  mutate(pc2=rollapply(plantc, 5, mean,fill=NA))

filr2 <- test2 %>% dplyr::filter(year == 2011) %>%
  dplyr::filter(ma2 > low_lai & ma2 < high_lai &
                  lai > low_lai & lai < high_lai)

filr2 <- inner_join(filr2, p, by='run')

these parameters were selected based on the sensitivity plots below.

plot time series of LAI

rather than aggregated, this plot is daily values. The red line is the mean, and blue lines are first and third quartiles.

sensitivity of 2011 LAI compared to 2017 LAI

I’ve been checking the LAI for these two years separately, to see if there is some parameter that matters more for the drought years leading up to 2017. Before your updates, epc.litter_moist_coef usually was high, not sure if that had to do with the litter and decomposition tht you changed. But these look pretty similar, which is a good thing!

compare range of RHESSys model output with mean values from data

finally, just plotting the distribution of the model output LAI in 2011 with the RS data. The model is underestimating LAI a bit. But overestimating plant carbon and height.

combine output and look at individual params

vis <- read.csv("/Volumes/GoogleDrive/My Drive/plant_params/DM_data/tree_and_turfgrass_monroe2_selected_weightedmeanVIs.csv")
quag_data <- vis %>% dplyr::filter(sp_code==
                                     "QUAG")
nd11 <- quag_data %>% dplyr::filter(year==2011) %>% dplyr::select(NDVI)
nd17 <-  quag_data %>% dplyr::filter(year==2017) %>% dplyr::select(NDVI)
ndvi_resil <- nd17/nd11

# sann <- sann %>% dplyr::filter(epc.root_growth_direction > 0.7 & epc.root_growth_direction <0.9)

ggplot(ndvi_resil) + geom_density(aes(x=NDVI)) +
  geom_density(data=sann, aes(x=resil, col='output'))

ggplot(sann) + geom_point(aes(x=epc.day_leafoff, y=resil))+geom_smooth(aes(x=epc.day_leafoff, y=resil, col=`2011`), method='lm') + ggtitle("day leafoff v. resil")
## `geom_smooth()` using formula 'y ~ x'

ggplot(sann) + geom_point(aes(x=sann$epc.day_leafoff, y=`2011`))+geom_smooth(aes(x=epc.day_leafoff, y=`2011`), method='lm') + ggtitle("day leafoff v. LAI in 2011")
## Warning: Use of `sann$epc.day_leafoff` is discouraged. Use `epc.day_leafoff`
## instead.
## `geom_smooth()` using formula 'y ~ x'

ggplot(sann) + geom_point(aes(x=`2011`, y=resil, col=epc.day_leafoff)) + ggtitle("LAI and resil") + geom_hline(aes(yintercept=mean(ndvi_resil$NDVI)))

ggplot(sann) + geom_point(aes(x=epc.min_leaf_carbon, y=resil, col=`2011`))+geom_smooth(aes(x=epc.min_leaf_carbon, y=resil), method='lm') + ggtitle("min_leaf_carbon v. resil")
## `geom_smooth()` using formula 'y ~ x'

ggplot(sann) + geom_point(aes(x=epc.storage_transfer_prop, y=resil, col=`2011`))+geom_smooth(aes(x=epc.storage_transfer_prop, y=resil), method='lm') + ggtitle("storage transfer prop v. resil")
## `geom_smooth()` using formula 'y ~ x'

ggplot(sann) + geom_point(aes(x=epc.resprout_leaf_carbon, y=resil, col=`2011`))+geom_smooth(aes(x=epc.resprout_leaf_carbon, y=resil), method='lm') + ggtitle("resprout leaf c v. resil")
## `geom_smooth()` using formula 'y ~ x'

ggplot(sann) + geom_point(aes(x=epc.root_distrib_parm, y=resil, col=`2011`))+geom_smooth(aes(x=epc.root_distrib_parm, y=resil), method='lm') + ggtitle("root distrib parm v. resil")
## `geom_smooth()` using formula 'y ~ x'

ggplot(sann) + geom_point(aes(x=`2011`, y=resil, col=epc.root_distrib_parm)) + ggtitle("LAI and resil") + geom_hline(aes(yintercept=mean(ndvi_resil$NDVI)))